小鼠

In [1]:
%matplotlib notebook
In [2]:
import numpy as np
import matplotlib as mpl
import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
import scipy.signal
import scipy.fftpack

mpl.rcParams["figure.dpi"] = 300
mpl.rcParams["figure.figsize"] = (10, 5)
In [31]:
import gc
gc.collect()
Out[31]:
0
In [35]:
waveform0 = plt.imread("mouse1-processed.png")
waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
waveform
Out[35]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [40]:
np.max(waveform), np.min(waveform)
Out[40]:
(0.99999999999999989, 0.0)
clear all
clc
tic
%图片读取
waveform0=rgb2gray(imread('Image1.bmp'));
waveform0=im2double(waveform0);
waveform0=waveform0(1:645,12:1024);
In [29]:
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
Out[29]:
<matplotlib.colorbar.Colorbar at 0x7f2230bf0550>
In [69]:
waveform.shape
Out[69]:
(644, 1012)

图片尺寸是644x1012

In [26]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
Out[26]:
<matplotlib.colorbar.Colorbar at 0x7f2223e0bc50>
In [10]:
fs = 300e+6
c0 = 1500
x0 = c0 / fs
t0 = 1 / fs
n = 30
nn = 2 * n + 1

fs=300e6;              %采样频率200MHz  中心频率80MHz
c0=1500;               %组织声速
x0=c0/fs;              %每个采样点的距离
t0=1/fs;               %采样间隔
n=30;                  %采样点数目2n+1
nn=2*n+1;
In [ ]:
slope = np.zeros((waveform.shape[0] - 2 * n, waveform.shape[1]))
intercept = np.zeros((waveform.shape[0] - 2 * n, waveform.shape[1]))
f = list(range(0, n+1)) # f = [0, 1, 2, ..., n]
for i in range(0, waveform.shape[1] - 2):
    print(i)
    for j in range(0, waveform.shape[0] - nn):
        temp = waveform[j: j + 2 * n + 1, i]
        temp = temp * np.hamming(nn)
        fft_temp = np.abs(np.fft.fft(temp))
        FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#         FT = fft_temp / np.max(fft_temp)
        FT1 = FT[0: n + 1]
#         print(FT1)
        RE_a, RE_b = np.polyfit(f, FT1, 1)
        slope[j, i] = RE_a
        intercept[j, i] = RE_b

slope, intercept

slope=zeros(RF_Row-2*n,RF_Column);
intercept=zeros(RF_Row-2*n,RF_Column);
f=1:n+1;
for i=1:RF_Column-2
    for j=1:RF_Row-nn
        temp=waveform(j:j+2*n,i);
        temp=temp.*hamming(nn);
        fft_temp=abs(fft(temp));
        FT=20*log10(fft_temp./max(max(fft_temp)));
        FT1=FT(1:n+1);
        [RE_a,RE_b]=line_reg1(FT1,f);
        slope(j,i)=RE_a;
        intercept(j,i)=RE_b;
    end
end
toc

代码优化

In [11]:
f = list(range(0, n+1))
def slopeFilter(X):
#     temp = X * np.hamming(nn)
#     fft_temp = np.abs(np.fft.fft(temp))
    fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
#     fft_temp = np.abs(scipy.fftpack.fft(X * np.hamming(nn)))
    FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#     print(X)
    return np.polyfit(f, FT[0: n+1], 1)[0]

# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2
In [21]:
%%timeit
slopeFilter(waveform[0: nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in true_divide
  

401 µs ± 26 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

slopeFilter需要约$400\,{\rm \mu s}$做一个窗口的计算。图片的大小是644x1012。所以总共有 $$ (644 - 61 + 1) \times 1012 = 591008 $$ 个窗口。预计整张图片需要 $$ 591008 \times 400 \cdot 10^{-6} = 236.4032 \,{\rm s} $$ 才能处理完。

希望能把单个窗口的运算时间压缩到$10 \,{\rm \mu s}$。这样每张图片需要的时间是 $$ 591008 \times 10 \cdot 10^{-6} = 5.91008 \,{\rm s} $$

In [29]:
%%timeit
for i in range(0, waveform.shape[0] - nn): 
    slopeFilter(waveform[i: i + nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in true_divide
  

220 ms ± 22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

每列需要处理$644 - 61 + 1 = 584$个窗口。按照理论需要 $$ 584 \times 400 \cdot 10^{-6} = 0.2336 \,{\rm s} $$ 才能结束。然而实际测得需要$200 \,{\rm ms}$。

如果按照实际测得的时间算的话，整张图片需要 $$ 1012 \times 200 \cdot 10^{-3} = 202.4 \,{\rm s} $$ 来处理。

In [9]:
%%time
slope2 = np.zeros((waveform.shape[0] - nn + 1, waveform.shape[1]))
f = list(range(0, n+1))
for i in range(0, waveform.shape[1]): # column
#     print(i)
    for j in range(0, waveform.shape[0] - nn): # row
        slope2[j, i] = slopeFilter(waveform[j: j + nn, i])
# slopeFilter(waveform[0: nn, 0])
slope2
# %time

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  

CPU times: user 2min 49s, sys: 38.3 ms, total: 2min 50s
Wall time: 2min 50s

用时$227\,{\rm s}$。

分析slopeFilter的延时

In [94]:
%prun slopeFilter(waveform[0: nn, 0])

 
/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  

173 function calls in 0.001 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.fft.fftpack_lite.cfftf}
        1    0.000    0.000    0.001    0.001 <ipython-input-11-84a0b36be0fd>:2(slopeFilter)
        2    0.000    0.000    0.000    0.000 {built-in method numpy.linalg.lapack_lite.dgelsd}
        1    0.000    0.000    0.000    0.000 linalg.py:1813(lstsq)
        1    0.000    0.000    0.000    0.000 polynomial.py:398(polyfit)
        1    0.000    0.000    0.000    0.000 function_base.py:3484(hamming)
       13    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.array}
        3    0.000    0.000    0.000    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        2    0.000    0.000    0.000    0.000 iostream.py:195(schedule)

发现np.fft比较慢，要$30 \,{\rm \mu s}$左右。如果用scipy.fftpack.fft代替的话，能快出$10 \,{\rm \mu s}$左右。

np.polyfit也是很慢的，要$180 \,{\rm \mu s}$左右。

替换掉numpy.fft改用scipy.fftpack.fft

In [67]:
f = list(range(0, n+1))
def slopeFilter(X):
#     temp = X * np.hamming(nn)
#     fft_temp = np.abs(np.fft.fft(temp))
#     fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
    fft_temp = np.abs(scipy.fftpack.fft(X * np.hamming(nn)))
    FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#     print(X)
    return np.polyfit(f, FT[0: n+1], 1)[0]

# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2

所以这里把np.fft.fft替换成scipy.fftpack.fft。

In [70]:
%%timeit
slopeFilter(waveform[0: nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in true_divide
  import sys

323 µs ± 7.95 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [73]:
%%timeit
for i in range(0, waveform.shape[0] - nn): 
    slopeFilter(waveform[i: i + nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in true_divide
  import sys

191 ms ± 3.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

果然快了一点。

事先生成好Hamming窗

In [74]:
%%timeit
np.hamming(61)

16.1 µs ± 786 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

每次生成Hamming窗都要花$16 \,{\rm \mu s}$。这部分时间可以省下来的。

In [5]:
f = list(range(0, n+1))
hammingWindow = np.hamming(nn)
def slopeFilter(X):
#     temp = X * np.hamming(nn)
#     fft_temp = np.abs(np.fft.fft(temp))
#     fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
    fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
    FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#     print(X)
    return np.polyfit(f, FT[0: n+1], 1)[0]

# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2

这里事先就生成一次Hamming窗，不再每次再在函数调用里生成一次。

In [120]:
%%timeit
slopeFilter(waveform[0: nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  

301 µs ± 7.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [81]:
%%timeit
for i in range(0, waveform.shape[0] - nn): 
    slopeFilter(waveform[i: i + nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  

170 ms ± 5.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

又省了$20 \,{\rm ms}$每列。

禁用多线程

In [91]:
np.ALLOW_THREADS = 0
In [97]:
%%timeit
slopeFilter(waveform[0: nn, 0])

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  

307 µs ± 16.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [ ]:
%%timeit
for i in range(0, waveform.shape[0] - nn): 
    slopeFilter(waveform[i: i + nn, 0])

没有什么作用。

利用多核

In [8]:
import ipyparallel as ipp
rc = ipp.Client()

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipyparallel/client/client.py:458: RuntimeWarning: 
            Controller appears to be listening on localhost, but not on this machine.
            If this is true, you should specify Client(...,sshserver='you@benjamin-Lenovo-IdeaPad-Yoga-13')
            or instruct your controller to listen on an external IP.
  RuntimeWarning)
In [9]:
lview = rc.load_balanced_view()
# result.display_outputs()
lview
Out[9]:
<LoadBalancedView None>
In [10]:
%%time
def task(i):
    i, waveform = i
    import numpy as np
    import scipy.fftpack
    
    fs = 300e+6
    c0 = 1500
    x0 = c0 / fs
    t0 = 1 / fs
    n = 30
    nn = 2 * n + 1
    
#     f = list(range(5, 29))
    f = list(range(0, n+1))
    hammingWindow = np.hamming(nn)
    slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    
    def slopeFilter(X):
        fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
        FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#         return np.polyfit(f, FT[5: 29], 1)
        return np.polyfit(f, FT[0: n+1], 1)
    
    for j in range(0, waveform.shape[0] - nn):
        slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
        
    return slope3, intercept3
    
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
#     results.append(lview.apply(task, i, waveform[:, i]))
#     print(i)

results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
    
results

CPU times: user 720 ms, sys: 62.8 ms, total: 783 ms
Wall time: 748 ms
In [7]:
results.serial_time
Out[7]:
448.07730399999934
In [8]:
results.wall_time
Out[8]:
116.608309

这是开4 workers的情况。用了$127 \,{\rm s}$。

In [97]:
results.serial_time
Out[97]:
322.10072999999977
In [98]:
results.wall_time
Out[98]:
166.383029

这是开2 workers的情况。用了$166 \,{\rm s}$。

In [9]:
slope3 = np.hstack([i[0] for i in results.result()])
slope3
Out[9]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [10]:
intercept3 = np.hstack([i[1] for i in results.result()])
intercept3
Out[10]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [13]:
%time scipy.ndimage.generic_filter(waveform[0: 100, 0: 20], slopeFilter, size=(61, 1))

/home/benjamin/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in true_divide
  

CPU times: user 724 ms, sys: 1e+03 ns, total: 724 ms
Wall time: 724 ms
Out[13]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [15]:
intercept
Out[15]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [11]:
slope3.tofile("mouse1-slope.numpy")
intercept3.tofile("mouse1-intercept.numpy")

图像

In [146]:
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
rows = np.array(range(slope2.shape[0]))
cols = np.array(range(slope2.shape[1]))
rows, cols = np.meshgrid(cols, rows)
ax.plot_surface(rows, cols, slope2, cmap="jet")
ax.set_aspect("equal")
# fig.colorbar(slope2)
In [ ]:
slopeImage = mpl.imsh
In [34]:
plt.imshow(slope3, cmap="jet")
plt.colorbar()
# plt.show()
Out[34]:
<matplotlib.colorbar.Colorbar at 0x7f2222b30ef0>
In [35]:
plt.imshow(20 * np.abs(slope3), cmap="jet")
plt.colorbar()
# plt.show()
Out[35]:
<matplotlib.colorbar.Colorbar at 0x7f2222a5acc0>
In [36]:
plt.imshow(intercept3, cmap="jet")
plt.colorbar()
# plt.show()
Out[36]:
<matplotlib.colorbar.Colorbar at 0x7f2222984c88>
In [21]:
plt.imshow(np.abs(intercept3), cmap="jet")
plt.colorbar()
# plt.show()
Out[21]:
<matplotlib.colorbar.Colorbar at 0x7fa454e706d8>
In [37]:
plt.imshow(slope3 * 16 + intercept3, cmap="jet")
plt.colorbar()
# plt.show()
Out[37]:
<matplotlib.colorbar.Colorbar at 0x7f2222931860>

选取肿瘤区和正常区

In [85]:
colNormal = 600
rowNormal = 350

colTumor = 380
rowTumor = 320

rows = 80
cols = 100

rowOffset = 30
In [43]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
Out[43]:
<matplotlib.patches.Rectangle at 0x7fc61ee9c668>
In [52]:
plt.imshow(slope3, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
Out[52]:
<matplotlib.patches.Rectangle at 0x7fc61787a860>
In [86]:
plt.subplot(1, 3, 1)
plt.hist(slope3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")

plt.subplot(1, 3, 2)
plt.hist(intercept3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.hist((slope3 * 16 + intercept3)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope3 * 16 + intercept3)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
Out[86]:
Text(0.5,1,'midband')
In [53]:
plt.subplot(1, 3, 1)
plt.boxplot([
    slope3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    slope3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")

plt.subplot(1, 3, 2)
plt.boxplot([
    intercept3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    intercept3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.boxplot([
    (slope3 * 16 + intercept3)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    (slope3 * 16 + intercept3)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")

# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()

小鼠2

In [4]:
import scipy.io
In [5]:
namespace = scipy.io.matlab.loadmat("mouse2.mat")
slope2 = namespace["slope"]
intercept2 = namespace["intercept"]
In [29]:
waveform0 = plt.imread("mouse2-processed.png")
waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
waveform
Out[29]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [30]:
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
Out[30]:
<matplotlib.colorbar.Colorbar at 0x7f2618ff4908>
In [32]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
Out[32]:
([], <a list of 0 Text yticklabel objects>)
In [35]:
%%time
def task(i):
    i, waveform = i
    import numpy as np
    import scipy.fftpack
    
    fs = 300e+6
    c0 = 1500
    x0 = c0 / fs
    t0 = 1 / fs
    n = 30
    nn = 2 * n + 1
    
#     f = list(range(5, 29))
    f = list(range(0, n+1))
    hammingWindow = np.hamming(nn)
    slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    
    def slopeFilter(X):
        fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
        FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
#         return np.polyfit(f, FT[5: 29], 1)
        return np.polyfit(f, FT[0: n+1], 1)
    
    for j in range(0, waveform.shape[0] - nn):
        slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
        
    return slope3, intercept3
    
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
#     results.append(lview.apply(task, i, waveform[:, i]))
#     print(i)

results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
    
results

CPU times: user 2.49 s, sys: 194 ms, total: 2.68 s
Wall time: 3.54 s
In [36]:
slope2 = np.hstack([i[0] for i in results.result()])
slope2
Out[36]:
array([[ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       ..., 
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        , -0.3075976 ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        , -0.30471012],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ]])
In [37]:
intercept2 = np.hstack([i[1] for i in results.result()])
intercept2
Out[37]:
array([[ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       ..., 
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  1.45241244],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  1.43369425],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ]])
In [40]:
plt.imshow(slope2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
Out[40]:
([], <a list of 0 Text yticklabel objects>)
In [41]:
plt.imshow(intercept2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
Out[41]:
([], <a list of 0 Text yticklabel objects>)
In [42]:
plt.imshow(slope2 * 16 + intercept2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
Out[42]:
([], <a list of 0 Text yticklabel objects>)
In [38]:
scipy.io.matlab.savemat("mouse2.mat", {
    "slope": slope2,
    "intercept": intercept2,
    "midband": slope2 * 16 + intercept2
})

选择肿瘤区和正常区

In [6]:
colNormal = 600
rowNormal = 200

colTumor = 620
rowTumor = 300

rows = 80
cols = 100

rowOffset = 30
In [45]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.xticks([])
plt.yticks([])
Out[45]:
([], <a list of 0 Text yticklabel objects>)
In [46]:
plt.imshow(slope2 * 16 + intercept2, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.xticks([])
plt.yticks([])
Out[46]:
([], <a list of 0 Text yticklabel objects>)
In [9]:
plt.scatter(slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), c=(0, 0, 0, 0.2), marker=".")
plt.scatter(slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), c=(1, 0, 0, 0.2), marker=".")
Out[9]:
<matplotlib.collections.PathCollection at 0x7f99a8532dd8>
In [84]:
plt.subplot(1, 3, 1)
plt.hist(slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")

plt.subplot(1, 3, 2)
plt.hist(intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.hist((slope2 * 16 + intercept2)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope2 * 16 + intercept2)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
Out[84]:
Text(0.5,1,'midband')
In [49]:
plt.subplot(1, 3, 1)
plt.boxplot([
    slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")

plt.subplot(1, 3, 2)
plt.boxplot([
    intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.boxplot([
    (slope2 * 16 + intercept2)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    (slope2 * 16 + intercept2)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")

# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()
In [47]:
plt.imshow((slope2 < -0.8) & (intercept2 < -11.25) & ((slope2 * 16 + intercept2) < -25), cmap="gray")
plt.colorbar()
plt.xticks([])
plt.yticks([])
Out[47]:
([], <a list of 0 Text yticklabel objects>)

血管

In [3]:
import scipy.io
In [4]:
namespace = scipy.io.matlab.loadmat("Ivus_data.mat")
namespace
Out[4]:
{'Ivus_data': array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        ..., 
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.]]),
 '__globals__': [],
 '__header__': b'MATLAB 5.0 MAT-file, Platform: PCWIN64, Created on: Tue Mar 20 20:13:56 2018',
 '__version__': '1.0'}
In [5]:
# waveform0 = namespace["Rabbit_aorta_80M_Sep8_pa_RF_01"]
waveform0 = namespace["Ivus_data"]
# waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
# waveform = np.nan_to_num(np.log(np.abs(waveform0)))
waveform = waveform0
waveform
Out[5]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [18]:
waveform.max(), waveform.min()
Out[18]:
(0.20000000000000001, -0.20000000000000001)
In [36]:
waveform.shape
Out[36]:
(3021, 3021)
In [1]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-ce43377d2b9b> in <module>()
----> 1 plt.imshow(waveform, cmap="jet")
      2 plt.colorbar()
      3 # plt.show()

NameError: name 'plt' is not defined
In [2]:
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-a370838f4a9c> in <module>()
----> 1 plt.imshow(waveform, cmap="gray")
      2 plt.colorbar()
      3 # plt.show()

NameError: name 'plt' is not defined
In [9]:
%%time
def task(i):
    i, waveform = i
    import numpy as np
    import scipy.fftpack
    
    fs = 300e+6
    c0 = 1500
    x0 = c0 / fs
    t0 = 1 / fs
    n = 30
    nn = 2 * n + 1
    
    f = list(range(5, 29))
    hammingWindow = np.hamming(nn)
    slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
    
    def slopeFilter(X):
        fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
        FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
        return np.polyfit(f, FT[5: 29], 1)
    
    for j in range(0, waveform.shape[0] - nn):
        slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
        
    return slope3, intercept3
    
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
#     results.append(lview.apply(task, i, waveform[:, i]))
#     print(i)

results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
    
results

CPU times: user 6.26 s, sys: 487 ms, total: 6.74 s
Wall time: 7.7 s
In [10]:
slope4 = np.hstack([i[0] for i in results.result()])
slope4
Out[10]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [11]:
intercept4 = np.hstack([i[1] for i in results.result()])
intercept4
Out[11]:
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
In [3]:
plt.imshow(slope4, cmap="jet")
plt.colorbar()
# plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-3-b3b93e8ecb1a> in <module>()
----> 1 plt.imshow(slope4, cmap="jet")
      2 plt.colorbar()
      3 # plt.show()

NameError: name 'plt' is not defined
In [4]:
plt.imshow(intercept4, cmap="jet")
plt.colorbar()
# plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-831ce4d62add> in <module>()
----> 1 plt.imshow(intercept4, cmap="jet")
      2 plt.colorbar()
      3 # plt.show()

NameError: name 'plt' is not defined
In [5]:
plt.imshow(slope4 * 16 + intercept4, cmap="jet")
plt.colorbar()
# plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-3592e6254b48> in <module>()
----> 1 plt.imshow(slope4 * 16 + intercept4, cmap="jet")
      2 plt.colorbar()
      3 # plt.show()

NameError: name 'plt' is not defined
In [32]:
scipy.io.matlab.savemat("rabbit.mat", {
    "slope": slope4,
    "intercept": intercept4,
    "midband": slope4 * 16 + intercept4
})

兔

In [3]:
import scipy.io
In [31]:
namespace = scipy.io.matlab.loadmat("rabbit_aorta_80MHz.mat")
namespace
Out[31]:
{'__globals__': [],
 '__header__': b'MATLAB 5.0 MAT-file, Platform: PCWIN64, Created on: Wed Apr 18 21:14:15 2018',
 '__version__': '1.0',
 'image_log': array([[-17.68643622, -28.67609654, -22.32602774, ..., -27.70455119,
         -20.26901032, -22.21793313],
        [-20.39598036, -20.73648393, -20.93056663, ..., -18.14376519,
         -19.28679759, -16.50690175],
        [-10.30797713, -11.67582775, -10.72473036, ...,  -8.56597864,
         -10.55324122, -10.2685423 ],
        ...,
        [-11.08561572,  -9.245724  , -12.79251056, ..., -13.47482297,
         -13.59563611, -10.87131917],
        [ -8.12298501, -23.97745454, -39.23335071, ..., -13.96027262,
         -15.47194327, -10.79594562],
        [-15.17570577, -17.85875206, -18.08302249, ..., -15.12017247,
         -17.02972531,  -9.71389133]])}
In [32]:
# waveform0 = namespace["Rabbit_aorta_80M_Sep8_pa_RF_01"]
waveform0 = namespace["image_log"]
# waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
# waveform = np.nan_to_num(np.log(np.abs(waveform0)))
waveform = waveform0
waveform
Out[32]:
array([[-17.68643622, -28.67609654, -22.32602774, ..., -27.70455119,
        -20.26901032, -22.21793313],
       [-20.39598036, -20.73648393, -20.93056663, ..., -18.14376519,
        -19.28679759, -16.50690175],
       [-10.30797713, -11.67582775, -10.72473036, ...,  -8.56597864,
        -10.55324122, -10.2685423 ],
       ...,
       [-11.08561572,  -9.245724  , -12.79251056, ..., -13.47482297,
        -13.59563611, -10.87131917],
       [ -8.12298501, -23.97745454, -39.23335071, ..., -13.96027262,
        -15.47194327, -10.79594562],
       [-15.17570577, -17.85875206, -18.08302249, ..., -15.12017247,
        -17.02972531,  -9.71389133]])
In [6]:
waveform.max(), waveform.min()
Out[6]:
(36.05074837859364, -62.46266160516997)
In [25]:
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
Out[25]:
<matplotlib.colorbar.Colorbar at 0x7f883a01e400>
In [13]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
Out[13]:
<matplotlib.colorbar.Colorbar at 0x7f8840121518>
In [13]:
slope5 = np.fromfile("rabbit-slope.numpy").reshape(-1, 1000)
intercept5 = np.fromfile("rabbit-intercept.numpy").reshape(-1, 1000)
midband5 = slope5 * 16 + intercept5
In [11]:
slope5 = np.hstack([i[0] for i in results.result()])
slope5
Out[11]:
array([[-0.8281674 , -0.83519525, -1.05915205, ..., -0.74537372,
        -0.75297503, -0.92240071],
       [-0.79317912, -0.82555709, -1.04656068, ..., -0.69080656,
        -0.71536978, -0.93683651],
       [-0.77144668, -0.81739586, -0.9697593 , ..., -0.61810374,
        -0.6847368 , -0.9681786 ],
       ...,
       [-0.65338191, -0.53003156, -1.05368961, ..., -0.56395875,
        -0.583135  , -0.8216247 ],
       [-0.67377837, -0.53368256, -1.00574435, ..., -0.56394465,
        -0.57963811, -0.82356964],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ]])
In [13]:
slope5.tofile("rabbit-slope.numpy")
In [12]:
intercept5 = np.hstack([i[1] for i in results.result()])
intercept5
Out[12]:
array([[-12.83554491, -11.27739111,  -9.73721505, ..., -15.77141538,
        -18.32175809, -15.7249512 ],
       [-12.63938482, -10.75036011,  -9.73047907, ..., -15.60803094,
        -18.30434427, -15.69724014],
       [-12.47890637, -10.29166449,  -9.81114317, ..., -15.99934698,
        -17.97963193, -15.63407603],
       ...,
       [-15.43360958, -18.00147884, -12.76840861, ..., -16.17695036,
        -17.9935765 , -16.93215717],
       [-15.37943776, -17.95404985, -13.02376167, ..., -16.24107044,
        -18.18744762, -16.79368366],
       [  0.        ,   0.        ,   0.        , ...,   0.        ,
          0.        ,   0.        ]])
In [14]:
intercept5.tofile("rabbit-intercept.numpy")
In [26]:
plt.imshow(slope5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
Out[26]:
<matplotlib.colorbar.Colorbar at 0x7f883a0a4518>
In [27]:
plt.imshow(intercept5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
Out[27]:
<matplotlib.colorbar.Colorbar at 0x7f8839e41f28>
In [28]:
plt.imshow(slope5 * 16 + intercept5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
Out[28]:
<matplotlib.colorbar.Colorbar at 0x7f883a14b780>

选择肿瘤区和正常区

In [7]:
colNormal = 100
rowNormal = 500

colTumor = 100
rowTumor = 200

rows = 80
cols = 100

rowOffset = 30
In [33]:
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
Out[33]:
<matplotlib.patches.Rectangle at 0x7f8839ca6160>
In [35]:
plt.imshow(slope5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
Out[35]:
<matplotlib.patches.Rectangle at 0x7f8839b6a0f0>
In [36]:
plt.imshow(intercept5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
Out[36]:
<matplotlib.patches.Rectangle at 0x7f8839a80ef0>
In [34]:
plt.imshow(slope5 * 16 + intercept5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
    plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
    plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
Out[34]:
<matplotlib.patches.Rectangle at 0x7f8839bc1f98>
In [25]:
plt.scatter(slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), c=(0, 0, 0, 0.2), marker=".")
plt.scatter(slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), c=(1, 0, 0, 0.2), marker=".")
plt.xlabel("slope")
plt.ylabel("intercept")
Out[25]:
Text(0,0.5,'intercept')
In [37]:
plt.subplot(1, 3, 1)
plt.hist(slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")

plt.subplot(1, 3, 2)
plt.hist(intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.hist((slope5 * 16 + intercept5)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope5 * 16 + intercept5)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
Out[37]:
Text(0.5,1,'midband')
In [38]:
plt.subplot(1, 3, 1)
plt.boxplot([
    slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")

plt.subplot(1, 3, 2)
plt.boxplot([
    intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")

plt.subplot(1, 3, 3)
plt.boxplot([
    (slope5 * 16 + intercept5)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
    (slope5 * 16 + intercept5)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")

# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()